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Abstract 

Turbulent backward-facing step flow was examined 
using four low turbulent Reynolds number k-e models and 
one standard high Reynolds number technique. A tunnel 
configuration of 1:9 (step height: exit tunnel height) was 
used. The models tested include: the original Jones and 
Launder, Chien; Launder and Sharma; and the recent 
Shih and Lumley formulation. The experimental refer- 
ence of Driver and Seegmiller was used to make detailed 
comparisons between reattachment length, velocity, pres- 
sure, turbulent kinetic energy, Reynolds shear stress, and 
skin friction predictions. The results indicated that the use 
of a wall function for the standard k-e technique did not 
reduce the calculation accuracy for this separated flow 
when compared to the low turbulent Reynolds number 
techniques. 

Introduction 

For the past two decades, researchers interested in 
evaluating wall bounded flows have relied upon two equa- 
tion turbulence models to facilitate this effort In the con- 
text of the Reynolds Averaged Navier Stokes (RANS) 
equations, this involved an assumption that the closure 
problem has been bridged by the Boussinesq approxi- 
mation. This approximation relates the Reynolds stress 
tensor to the velocity gradient tensor via an eddy vis- 
cosity. The family of two equation models defined this 
eddy viscosity in terms of characteristic turbulent veloc- 
ity and length scales which resulted from the solution of 
two additional partial differential equations. One example 
of this theory is the energy-dissipation rate (k-e) model. 
However, for wall bounded flows, the two turbulent trans- 
port equations behave differently for each separate region 
of near-wall flow. Modelling this near-wall flow behavior 
is crucial to the success of an internal flow simulation. 

One of the more challenging internal flows used to 
benchmark a given model is the backward-facing step 
(BFS). The presence of a strong recirculation bubble, a 
reattaching shear layer, a redeveloping channel flow and 
a relatively simple geometry have made this a particularly 
attractive flow to analyze; it has become a de facto 
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“complex canonical” flow 1 . BFS flow is an excellent 
candidate for examining the relationship of near-wall 
treatment to the k-e model performance for separated 
flows. 

Two definite schools of thought have emerged on 
how to treat this region of low turbulent Reynolds number 
(Ret) flow near a solid surface. The standard, or high 
Ret (HR) formulations utilized wall functions to model 
the near-wall flow, while the low Ret (LR) formulations 
incorporated the use of damping functions to model the 
effects of this region. 

At a first glance, the empiricism involved in the use 
of wall functions tends to make the more sophisticated 
LR models appear more attractive. Wall functions were 
based solely upon the behavior of attached equilibrium 
turbulent flow. There is, however, a significant cost 
savings to be realized from the appropriate use of HR 
models in terms of grid size. Alternatively, the LR 
models offer the ability to solve the transport equations 
through the viscous sublayer. Thus a straightforward 
boundary condition procedure is applied at the wall. This 
permits the near-wall velocity gradient to be found as 
a result of the calculation, not merely as an a priori 
assumption. It is this feature which has attracted such a 
great deal of attention to these LR models in the past two 
decades. However, the development of these models was 
also primarily based upon the characteristics of attached 
bo indary layer flow as well. 

Although much effort has recently been expended 
on the study of BFS flow with the k-e model, the over- 
whelming majority of analysis has utilized the wall func- 
tion approach due to the reduced numerical complexity. 
To the author’s knowledge, only So et al. 2 and Avva et 
al. 3 have examined the LR k-e formulation for BFS, and 
this was limited to the Chien model. While this model is 
known to be quite robust for a wide variety of flows, an 
anomalous behavior exists near the separation and reat- 
tachment points. This will be discussed in detail later. 
Several LR k-e models exist which do not demonstrate 
this same behavior for separated flow. 

The intention of this study is to investigate several 
LR models and compare against the performance of the 
HR formulation for separating flow. The first section 
describes the governing equations and turbulence models, 
as well as the experimental study of backward-facing step 
flow by Driver and Seegmiller 4 . The following section 
details the numerical methods used for this simulation. 
The final two sections include a discussion of the results 
and the conclusions of this research. 
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Background 


is expressed below. 


Governing Equations 

The RANS equations for incompressible flow are 
shown below in equation system 1. The dependent 
variables include the pressure and velocity as well as 
the Reynolds stress tensor, (uiuj). Here the quantities 
(P,U,V) denote the mean values while (p,u,v) are the 
fluctuating components. 


Ui,i = 0 

Ui.t + Uj Uij + ip.i - z/Ui.jj + (uluj) ^ = 0 

The Boussinesq approximation closes the system by as- 
suming the relation given in equations 2 . 


U,Uj 


-"«(Ui,j + Uj.i) + -k6ij 


k = 


1 

2 UlU| 


( 2 ) 


Notice that this is a function of the mean velocity gradi- 
ents and the turbulent kinetic energy, k. Now the eddy 
viscosity must be modelled. Gas kinetic theory yields a 
definition of the kinematic viscosity as a function of char- 
acteristic velocity and length scales. An eddy viscosity 
can be formulated if characteristic turbulent velocity and 
length scales can be postulated. For the k-e models, the 
velocity scale is given as the square root of k , while the 
length scale is a ratio involving k and the dissipation rate 
of turbulent kinetic energy, e. Thus, the isotropic turbu- 
lent viscosity is defined as v t = 

Transport equations are used to resolve the values of 
k anc e and are shown here in equation 3. 


Art+Uj k t j — 


v+ 


Ck 


= P-e + D 






=CxJ^V-C^+E 


( 3 ) 


where production of turbulent kinetic energy is defined 
as V- -(uTuJ)^(Uij + Uj,i). The constants C^, Ci, C 2 , 
<r k and ere are defined a priori. The terms D, E, f^, fi 
and f 2 are correction terms for the LR formulation and 
are necessary to give asymptotically correct behavior if 
the numerical domain is extended to the wall. 


k-e Model: HR Formulation 

In order to bridge the viscous sublayer, a two layer 
“law of the wall” was used to establish the values of kbc 
and £bc at the first point off the wall. The two layer model 
used in this test is similar to that given in reference 5 . For 
two dimensional boundary layer flow, the wall function 


„+ _ / y + for y + < 11 \ 

\ 5.5+ 2.5/n(y + ) for y + > 11 J 



Distance to the nearest wall was specified as (y) and the 
shear stress was denoted (r). If one assumed that this first 
interior grid point was properly positioned at the edge of 
the buffer layer, then empirical evidence suggests that the 
production of k was equal to the dissipation of kr. 


Pbc — £bc ( 5 ) 

The numerical boundary conditions for k^ and ebc were 
then derived from this equilibrium expression combined 
with the definition of u t and equations 4. 

% = th C = 2.5C4 — ( 6 ) 

U? y 

While this assumption is not valid for separated 
flow due the fact that the logarithmic overlap region 
is not present for this flow regime, we shall see later 
that the predictive capability of the HR formulation is 
not significantly diminished. To complete the model 
specification, the values of the constants and functions 
are listed in the tables below. 

k-e Model: LR Formulation 

Four different LR turbulence models were investi- 
gated in this study including: the original Jones and Laun- 
der (JL) 6 , Chien (CH) 7 , Launder and Sharma (LS ) 8 and 
the recent Shih and Lumley (SL ) 9 formulation. These 
models can be grouped into two broad categories accord- 
ing to the form of the near-wall corrections. The JL and 
LS models were developed as functions of the dependent 
variables alone, while the CH and SL models were based 
upon both the dependent and independent variables. All 
four models can be implemented in a domain extending 
to the wall. These are just a few of the many different 
models available, and the interested reader is referred to 
references 1<U1 * 12 for a more thorough review. 

The influence of the viscous dominated near-wall 
region was affected through the use of explicit damping 
functions and additional destruction terms. As mentioned 
earlier, these models were also developed to match the 
observed behavior of an attached boundary layer flow. 
Thus these LR formulations also suffer from an inherently 
empirical nature. The values of these functions, destruc- 
tion terms, and additional constants are listed in the tables 
below. Note that for the SL model, the damping function 

was derived in terms of Re k = ^ 7 ^, and the following 
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constants: ai = —1.5 x 10 -4 , a 3 = —1.0 x 10“ 9 , and 
a 5 = —5.0 x 10" 10 . The boundary conditions for the JL, 
LS, and CH models were specified as 

Cbc — kbc — 0 (7) 


while for the SL model the boundary conditions were 
slightly more complicated. 



Ci 

c 2 

Cm 

<r k 


JL 

1.55 

2.0 

0.09 

1.0 

1.3 

LS 

1.44 

1.92 

0.09 

1.0 

1.3 

CH 

1.35 

1.8 

0.09 

1.0 

1.3 

SL 

1.44 

1.92 

0.09 

1.0 

1.3 

HR 

1.44 

1.92 

0.09 

1.0 

1.3 


Table 1 The k-e model constants 




h 

<V 

JL 

1.0 

l-.3e-*‘? 

e ((H-R« t /50)) 

LS 

1.0 

l-.3e-*‘? 

e ((l + Re t /50)2 ) 

CH 

1.0 

Re ? 

1— .22e -_ 36" 

1 _ e -0.0115t,+ 

SL 

1.0 

Rc t 

1— ,22e s. - 

|"l g(aiRk+a3R^-|-a5n-k)| 2 

HR 

1.0 

1.0 

1.0 


Table 2 The LR k-e damping functions 



D 

E 

JL 

1 

to 

to 

2i/i/*UijkUijk 

LS 


2^(Uij k Ui >jk 

CH 

2 uk 

y 2 

t/ 3 e 

SL 

0 

^UijkUijk 

HR 

0 

0 


Table 3 The LR k-e additional terms 


Shih and Lumley have derived boundary conditions 
based upon an asymptotic behavior in the very near-wall 
region. As one approaches the wall from outside the 
boundary layer, a turbulent limit point is reached beyond 
which self sustaining turbulence cannot exist If this point 
is designated by the subscript a, then Shih and Lumley 
proposed the relations of equation 8. 

€« = 0.251^ ,ifc o = 0.25u? (8) 


For practical purposes, this limit point was assumed to 
extend to the wall; hence e bc — £ Q , he = k a . 

A comment on implementation in two dimensional 
flow is necessary for several models. The implementation 
of HR and SL involved the use of a value (y) which was 
based upon the minimum distance to a no-slip boundary 
in the numerical domain. This appeared to work quite 
well. For the CH model, the minimum y + value was used. 
This straight forward extension of the CH formulation to 
two dimensions resulted in a grid dependent solution in 
the region surrounding reattachment. 

BFS Flow: Experimental Configuration 

Many different experimental data exist for valida- 
tion of turbulent flow over a backward-facing step. The 
Stanford conference of 198Q-81 used the Kim et al. 13 
data for a 1:3 (step height: exit tunnel height) tunnel 
to benchmark model performance. The velocity profiles 
were measured via hot wire anemometry. However, tur- 
bulence data was unavailable in the recirculation zone. 
Recently, Kjelgaard 14 has examined backward-facing step 
flow using a three component laser velocimeter in a 1:2 
tunnel configuration. Mean velocity profiles and turbu- 
lent normal and shear stresses were tabulated throughout 
the flow domain. However the experiment of Driver and 
Seegmillei 4 has been chosen here because the measure- 
ments included skin friction, pressure, turbulent normal 
and shear stress, and velocity profiles from the step down- 
stream through the recovery region. A two component 
laser velocimeter was used to obtain the velocity and tur- 
bulence data. The tunnel expansion ratio of 1:9 (zero 
deflection angle) produced a much lower pressure gra- 
dient due to freestream sudden expansion than the other 
aforementioned tunnel geometries. 




Uref 



Fig. 1 Tunnel configuration for the 
backward-facing step of Driver and Seegmiller 4 

The tunnel configuration can be seen here in figure 
1. The inlet was 80H upstream from the expansion and 
the exit was downstream an additional 60H. The experi- 
mental Mach number was 0.128. The turbulent boundary 
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layer thickness was recorded as 1.5H at a location 4 step 
heights upstream of the expansion. The experiment was 
conducted at a Reynolds number (Ren) of 33420, based 
upon the step height According to Driver and Seeg- 
miller, this value insured that the boundary layer was 
fully turbulent as it passed over the step. 


Numerical Method 


RANS Equations 

A modified version of the finite volume code, 
DTNS2D 15 was chosen for this study. Chorin’s 16 pseudo 
compressibility technique was used to resolve the in- 
compressible RANS equations. The system of equations 
solved for the two dimensional pseudo compressibility 
method differs from the incompressible flow equations 
by the addition of a time dependent term in the conti- 
nuity equation, This is shown in nondimensional 

form for 2D laminar flow in equations 9. Note that x 
and y are the independent variables and Re refers to the 
Reynolds number. The variable (P) is defined as the pres- 
sure normalized by the constant density. For turbulent 
flow, the kinematic viscosity is replaced with the sum of 
turbulent and kinematic viscosities. This system can be 
regarded as hyperbolic in nature while the incompress- 
ible flow equations are elliptic. The pseudo sound speed, 
c = y/Uf 4- /?, is governed by the value of the parameter 
/? and the contravariant velocity U Cy whereas the physical 
sound speed is infinite. For this study, the value of /? 
was set to unity. 


b (q) + £ (fi+gi) + |^ (f2+sa) - 0 


■p" 


U/J 


V/? 

u 

, fl = 

U 2 + P 

, u = 

uv 

V 


uv 


V 2 + P 


-1 

" 0 " 

-1 

0 1 

au 

0U 

gl " Re 

j£ 

’ 62 “ Re 

ii 


- dx - 


L J 


The chief advantage of the pseudo compressibility 
technique is that the system of equations can now be 
marched in time towards a steady solution. With this 
in mind, the approximate factorization scheme was cho- 
sen to resolve the governing equations 9. The system 
is expressed in generalized coordinates (f (x, y), rj(x, y )) 
below in the equation system 10. 


^( q ) + ^ ( f » + «0 + b ( f2+g2 ) " 0 


j = 


d(t,v) 


q = 


d(x,y)' 

f _ £rfl + £y*2 ^ + Vy^2 

f 1 " J ’ ^ " J 


gl = 


_ Crgl + £yg2 


g2 = 


*7xgl + *7yg2 


(10) 


This system can be factored as shown below in equations 
11 


V d- + 0+ d‘ 

-T- + "qtA + -T— A + 7tA v 

.At di d£ d£ 


Aq* = RHS 


f V d~ . d+ & 1 V . 

{s.W + ^ B " + ^ B °} A4 = 37 Aq ' 

( 1 !) 

where A and B refer to the approximate convective flux 
Jacobians and A v and B v are the viscous flux Jacobians. 
The cell volume and timestep are defined as V and At, 
respectively. The time integration is complete given that 


RHS = 




( 12 ) 


and by defining the update q n+1 = q“ + Aq. The no- 
tation (d + , d~ , &) refers to the forward, backward, and 
central difference operators, respectively. The approxi- 
mate convective flux Jacobian A ± = RA 1 L is defined 
in terms of eigenvalues and the left and right eigenvec- 
tors, shown here for completeness. 


A = 


U c + a 
0 
0 


0 

U c 

0 


0 

0 , 
U c -a 


± _ (A±JAJ) 
2 


(13) 


The contravariant velocity is defined in terms of the 
metrics which correspond to both generalized coordinate 
directions (£, 77 ), depending on which Jacobian is being 
formed. 

U c = n x U + Ti y V 

N = n 2 + n 2 (14) 

« € (f , ri) 

The LHS was discretized to first order accuracy; how- 
ever, the RHS was based upon the second order accu- 
rate upwind TVD scheme of Chakravarthy and Osher 17 . 
The approximate Riemann solver and variable extrapola- 
tion technique are described in Gorski 15 . The interested 
reader is encouraged to examine this reference for fur- 
ther details. It was not necessary to utilize the limiter 
function in order to get a smooth solution for the RANS 
equations. 


Right Eigenvectors 

0 

1 

'cT 

1 

0 

—a(U c + a) 

-nJIJ c -a) + m 

n y 

nJU c +a)—\]N 

-^y(U c -a) + VN 

n x 

n^J e +a)-VN 


J 


4 


Table 4 The right eigenvectors of the 
2D incompressible RANS equations 


Left Eigenvectors 


1 

nxCCArhO 

n v (UA^) 

2 a 5 

2 a*N 

2a 2 N 

Un t -Vn, 

VU,+n,0 

U Ur+n,0 

a 2 

a 2 

a 2 

1 

2a 2 

n x(U c —a) 
2a 2 N 

n v (U c —a) 

2a 2 N 


Table 5 The left eigenvectors of the 
2D incompressible RANS equations 


Turbulent Transport Equations 

The approximate factorization technique was also 
applied to the turbulent transport equations 3. For two 
dimensional flow, the system is expressed as follows: 


|r( qt ) + + ^ (fzt+g ^ - st 


q‘ = 


k 
e 

t z! 

gl Re 


f * ~ 
I M — 




UJb 

Uf 


f 1 — 
*2 — 


Vk 

Ve 

dk 


,g2 


(15) 


The convective flux Jacobian for system 16 is equal 
to the identity matrix multiplied by the contravariant 
velocity, U c . The eigenvalue matrix is equivalent to the 
convective flux Jacobian matrix. Thus the left and right 
eigenvector matrices are both equal to the identity matrix. 
This permits a straightforward application of the same 
techniques used upon the flow equations 11. 

Michelassi and Shih 12 proposed a linearization for 
the source term Jacobian H l based upon retaining the 
positive leading order terms. After considerable effort, a 
similar formulation was found to be the most effective. 
These terms are shown in table 6. All off diagonal terms 
were set to zero so as not to impede diagonal dominance. 
The coefficients (a^Qrj) were used to include the implicit 
treatment of the source terms in both directions as needed 
for elliptic flows. The convergence was best when both 
were set equal to one half. The TVD variable extrap- 
olation is used for resolving the right hand side of this 
system, also. The limiter function is necessary for this 
system to avoid the nonphysical overshoots introduced 
by the linearization of H*. 


S l = -i- 
Re 


V — cRe-}-D 
CifiiP-C 2 f2T Ae + E 


with i/jc = v + v € = v + *£ 


and the production term 


v = *{2[«)’+ (f ) 2 ] + (f + Ato 

transforming the equations into generalized coordinates, 
the system was approximately factored as: 

{l? + | A ‘ ++ + | A * = RHS ‘ 
{^ + | B ‘ ++ ^ B *‘ + |; B »-“' H ‘W =>»'*• 

5 (fi + £i) d (^2+62) 


RHS‘ = — 


di 


drj 


+ S* 


( 16 ) 


Evaluation of the LR source terms for the CH model 
was straightforward. However, the other three LR models 
required an evaluation of the volume integrals: 


J J (Ujjk) 2 dt> = (U Uk ) 2 V 


//H,] — HJv 


(17) 


The following second order accurate linearizations were 
implemented: 

k 2 


(Ui,jk) 2 = (Ui, jk )'+0(Z\x 3 ,Ay 3 ) 


(Vk) . = f [(VT) J J +0{Ax\ Ay 3 ) 


(18) 


model 

dS\ 

dS\ 


dk 

de 

JL 

t 

k 

-c 2 hi 

LS 

€ 

k 

-C*hi 

CH 

2 

E 


Be t y 2 

£ 

SL 

£ 

k 

-2 C 2 hi 

HR 

0 

0 


Table 6 Diagonal terms for the 
approximate source Jacobian, H l 

BFS Flow: Numerical Configuration 

The numerical domain extended 15H upstream and 
60H downstream of the step. These values were found 
to generate results independent of the inlet and exit lo- 
cations, respectively. A mesh is constructed from three 
individual blocks with shared boundaries. However, the 
grid dimensions will be referred to as the total number of 
cells in each direction over the entire domain. Recently, 
Thangam and Hur 18 found that for a 1:3 tunnel configura- 
tion, a 166x73 cell Cartesian grid was necessary to fully 
resolve the flow features using a HR model. This grid was 
stretched towards the step in the streamwise (x) direction 
and uniform in the transverse (y) direction. However, due 
to the different nature of this 1:9 tunnel configuration, a 
separate grid independence study was carried out for the 
HR formulation. 

Three grids were used in the mesh dependence study. 
Due to the assumptions of the wall function given in 
equation 4, the normal distance from the wall to the first 
cell center was specified to yield y+ « 30 for regions 
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outside the recirculation zones. The three grids exam- 
ined were 62x40, 122x80, and 250x160 cells. Figure 2 
demonstrates the grid independence of the medium sized 
mesh results using profiles of the velocity and turbulence 
quantities at five step heights downstream of the step. By 
examining just the velocity data, one might assume that 
grid independence is achieved using the coarsest grid. 
However, the turbulence profiles reveal an improvement 
between the coarse and medium mesh results. The dif- 
ferences between the medium and fine mesh predictions 
for all quantities are negligible. Thus, the medium grid 
of 122x80 cells was chosen to conduct the HR testing 
on. For the LR formulations, the near-wall grid spacing 
was dictated by a discretization criterion of r/f « 1.0 
for the first cell along a no-slip boundary. Furthermore, 
15-20 cells were placed below the y + « 30.0 location. 
The LR grid was adapted using Turbol 19 so that the thin 
layer of high aspect ratio cells, necessary to resolve the 
separating shear layer, did not extend through the whole 
numerical domain. The mesh configurations used can be 


seen in figure 3 below. 


° Driver and Seegmiller 

- 62x40 cell mesh 

- 122x80 cell mesh 

- 250x160 cell mesh 



Fig. 2 Mesh refinement study: 
velocity and turbulence profiles 
downstream of the step at x/H = 5. 


HR Grid 


22x58 Block 1 
100x58 Bloci2 
100x22 Block3 



LR Odd 49x64 Block 1 

169 x64 Block 2 

169x50 Block 3 



Fig. 3 Multiblock mesh configuration for the HR and LR k-t models 


6 


Experimental data for the mean velocity and turbu- 
lent intensity profiles are available for the inlet tunnel 
at a location four step heights upstream of separation. 
These values were used at the numerical inlet. Pitot tube 
measurements of the velocity profiles helped to insure 
a reasonable massflow rate at the inlet. This was criti- 
cal to proper resolution of the pressure drop across the 
tunnel. However, the boundary layer along the wall op- 
posite the step was not accurately resolved using the laser 
Doppler velocimeter. Thus, the profiles of the turbulent 
quantities were mirrored about the inlet tunnel centerline. 
Even though the conditions were applied at fifteen step 
heights upstream of separation, the effect upon the result 
is negligible. Figure 4 demonstrates that the comparison 
between experiment and calculation at four step heights 
upstream from separation is good. 



Fig. 4 Inlet tunnel conditions at x/H=-4.0 

The velocities were specified and pressure was ex- 
trapolated at inlet, and vice versa at exit For the tur- 
bulence equations, both k and c were specified at inlet 
and extrapolated at exit Using the two turbulent inten- 
sity profiles, a kinetic energy profile was estimated using 
k = f(uu + vv). This implied that the third component 
was an average of the other two. Kjelgaard’s results 14 
support this assumption. The value for the mixing length 
(lm) was calculated using the simple ramp function 20 de- 
scribed in figure 5. 



Fig. 5 Ramp function definition 
of turbulent length scale 


Results 


The primary reattachment length is often used for 
measuring the overall performance for backward-facing 
step flow calculations. The experimental reattachment 
length for this tunnel geometry was given by Driver 
and Seegmiller as 6.26H. Numerically, this value was 
extrapolated from the velocity field and corresponds to 
the location where the value of wall shear stress was 
zero. The results are shown here in table 7. The LS, 
CH and HR performed the best, followed by the SL and 
JL models. 


model 

JL 

LS 

CH 

SL 

HR 

x r 

4.9 

5.4 

5.4 

5.1 

5.5 


Table 7 Primary reattachment length 


The velocity fields predicted by all five models were 
nearly identical. The superimposed profiles of figure 6 
clearly indicate this. Except for the reattachment region, 
the velocity profiles were indistinguishable. They have 
been plotted merely to demonstrate this point The wall 
static pressure coefficient, shown in figure 7 below, was 
more helpful in discriminating between the various model 
results. The trends evident in the reattachment length data 
were echoed here for the pressure rise along the stepside 
wall, with the exception of the Chien result It should be 
noted that all five formulations underpredicted the pres- 
sure drop in the recirculating flow and overpiedicted the 
pressure rise near reattachment The CH model underpre- 
dicted the pressure rise more than any other technique, 
while the LS model yielded the best agreement with the 
experimental data. The HR model has an especially dif- 
ficult problem in predicting the proper coefficient in the 
redeveloping flow region, although all four LR models 
overpredict this value to some degree as well. This might 
be a result of a slight error in massflow rate. This offset 
is also visible in the numerical results of Sindir 4 . 
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o Driver and Seegmiller 

Jones and Launder 

Launder and Sharma 

Chien 

Shih and Lumley 

Standard High Re, 




Fig. 6 Mean velocity profiles downstream of the step 


o Dover and Seegmilter. stepside wall 
a Driver and Seegmiller. opposite wan 

Jones and Launder 

Launder and Sharma 

Chien 

Shih and Lumley 

Standard High Re, 



Fig. 7 Wall static pressure coefficient downstream of the step 
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The profiles for the turbulent kinetic energy are also 
shown below. In the neighborhood of reattachment, the 
peak values were consistently located too far from the 
wall. In figure 8 the profiles at = 1 agree in location 
and magnitude for the peak value of k. Behavior within 
the recirculation bubble was best predicted here by the 
Chien model. The data at ^ = 2.5 appeared the most 
irregular and it is difficult to say which model performs 
best. The predictions appear to vary the most at this 
location. The turbulent kinetic energy was not well 
predicted by any technique for the next three profiles 
shown. However, the differences between the five k-e 
models are becoming much smaller. The peak values 


were also slightly overestimated. The remaining five 
profiles, beginning with ~ = 10, describe the recovery 
of redeveloping channel flow. The models overestimated 
the turbulent kinetic energy of the reattached shear layer 
initially, but this was significantly improved after ^ = 
20. For the HR model, the peak values of k are slightly 
lower throughout the redeveloping flow region. It should 
be noted that overall, the standard HR k-e model predicts 
similar results to the LR formulations for the turbulent 
kinetic energy, except in the region of high gradient near 
the separation line. This can be most easily seen in the 
first four profiles near y/E = 1.0. 


o Driver and Seegmiller 

Jones and Launder 

Launder and Sharma 

Chien 

Shih and Lumley 

Standard High Re, 




k/(U 2 fe( x10 2 ) 


Fig. 8 Turbulent kinetic energy profiles downstream of the step 
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Similar behavior was observed in the prediction of 
the turbulent shear stress shown in figure 9. Relatively 
good agreement was exhibited for all models at the § = 1 
location. However, by -g = 2.5 the peak value of turbu- 
lent shear was overpredicted by the JL and CH models. 
Also, the vertical position was located too far from the 
wall. It is interesting to note that the various LR profiles 
are bracketed here by the maximum values of JL and 
the minimum of LS. These two differ only in the damp- 


ing function fy* and the constants Cj and C 2 . As one 
proceeded downstream towards reattachment, the differ- 
ences between models were reduced and the overpredic- 
tions more pronounced. After reattachment, the various 
k-e models continued to overpredict the turbulent shear 
stress even throughout the recovery zone, although this 
effect was reduced as the exit was approached. Note that 
the HR results were predicting peak values slightly lower 
than the LR models. 


o Driver and Seegmiller 

Jones and Launder 

Launder and Sharma 

Chien 

Shih and Lumley 

Standard High Re, 




Fig. 9 Turbulent shear stress profiles downstream of the step 


For the backward-facing step flow, prediction of the 
skin friction coefficient along the stepside wall was the 
most difficult task. Results for the five models tested 
are given below in figure 10. The HR formulation pre- 
dicted the wall shear stress better than any of the LR 
formulations tested. However, all of the models have in- 
accurately predicted the near-wall velocity gradient for 
the recirculating flow close to reattachmenL The Chien 
model yielded the largest overshoot in Cy. This was in 
contrast to the relatively good reattachment length predic- 


tion of table 7. However, an examination of the Chien 
damping function indicated a dependence upon the 
y + parameter, or the nearest wall shear stress. In the 
neighborhood of a separation or reattachment point, an 
anomalous region of very low eddy viscosity results due 
to excessive damping in the neighborhood of zero wall 
shear stress. Figure 1 1 below demonstrates the effect of 
this y + dependence on the calculation of eddy viscosity. 
The behavior of the SL model in figure 12 was typical 
of the JL, LS, and HR results. 
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o Driver and Seegmilier 

Jones and Launder 

Launder and Sharma 

Chien 

Shih and Lumley 

Standard High Re. 



Fig. 10 Friction coefficient for the stepside wall downstream of the step 



Fig. 11 Nondimensional eddy viscosity field for the CH model 



Fig. 12 Nondimensional eddy viscosity field for the SL model 
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Close examination of the experimental skin friction 
data indicates a small comer eddy was present at ~ = 0.5. 
However, by ^ = 1.8 the primary recirculation zone was 
encountered. Numerical evidence of the comer eddy was 
demonstrated by two of the LR models: CH and SL. The 
JL, LS and HR formulations did not reveal any secondary 
recirculation zones. See table 8 below. This zone was 
confined to a very small region in the comer where the 
Ret values are low. 

A trend in the length of this secondary recirculation 
zone can be correlated to the damping function used in the 
definition of eddy viscosity for the LR models. Following 
the example of Patel et al. 11 , the damping functions 
were examined for attached boundary layer flows. From 
figure 13 we can see that for the JL and LS models, 
the damping in the viscous sublayer was never complete 
because the functions tend to small positive values instead 
of zero. In addition, there is a steep rise in the damping 
function around % 10. The SL and CH models both 
damped the eddy viscosity to zero at the wall, and also 
diffused the viscous effects further out from the wall. The 
most diffusive of the near- wall treatments was indeed the 
CH model. The largest secondary recirculation bubble 
also corresponds to the CH calculation. It appears that 
the interplay between the near-wall, viscous dominated 
flow and the fully turbulent, high Re; flow is critical to 
proper resolution of this comer region. However, more 
experimental and numerical work needs to be undertaken 
before any meaningful conclusions can be drawn. 


model 

JL 

LS 

CH 

SL 

Zr2 

0.0 

0.0 

0.8 

0.1 


Table 8 Secondary flow reattachment length 


Jones and L^undor 

— Launder and Sftarma 

Chwn 

SNh and Lurrtay 



y 


Fig. 13 Damping function behavior 
for attached boundary layer flow 

The preceding calculations were carried out on the 
Cray YMP at NASA Lewis. On average, a LR calculation 


needed roughly 5000 iterations to converge (CFL=3.0), 
depending upon the quality of the initial condition. Like- 
wise, a HR calculation required approximately 3300 iter- 
ations to converge to the same extent The convergence 
criterion was based upon the L 2 norm of the transient 
terms dropping five orders in magnitude. This was real- 
ized in all but the LS model, which tended to level off 
after a drop of three and a half orders. 


Conclusion 

Results for flow over the backward-facing step tun- 
nel configuration of Driver and Seegmiller were calcu- 
lated for the standard HR, and four LR Jc-e models. The 
velocity profiles and pressure coefficient data demon- 
strated how similar the results were among all five for- 
mulations. Slight differences appeared in the vicinity of 
reattachment All five models underpredicted the rede- 
veloping flow velocity profiles. Similarly, all five models 
yielded an early prediction of the pressure rise. A more 
discriminating measure of the near- wall velocity field was 
the reattachment length, x r . This was tabulated in table 
7. It appeared that the LS, CH, and HR models yielded 
similar performances in underpredicting this parameter 
by approximately 14%. Examination of the experimental 
velocity data near reattachment indicates that the reversed 
flow is confined to a very thin near-wall region; perhaps 
an improved damping strategy can be devised to resolve 
this low Ret region. However, I believe that the accurate 
prediction of skin friction is the most significant challenge 
to the LR k-e models for backward-facing step flow. 

Only the HR technique can claim to predict the wall 
shear stress in a bounded sense. The SL model produced 
the best prediction of the LR techniques, while the CH 
model produced the least reliable wall shear stress pre- 
diction. This appeared to be a result of the y * based 
damping function, I believe that a modification to 
the CH damping function which reduces this local sep- 
aration/reattachment effect may improve the shear stress 
prediction. 

The turbulent kinetic energy and shear stress profiles 
were generally overpredicted in the recirculation zone. 
The position of the peak values were consistently lo- 
cated too far from the wall as well. Perhaps the strong 
anisotropy of the flow contributes to this. Two anisotropic 
eddy viscosity models 5,21 have been reported to improve 
behavior for this flow. 

The significant conclusion of this study is that the 
assumption of local equilibrium inherent in the wall func- 
tion technique yields a credible result for separated flow. 
This conclusion is shared by So et al. 2 Furthermore, the 
near- wall modelling techniques employed in the JL, LS, 
SL, and CH appear to be much more sensitive to adverse 
pressure gradients than does the HR formulation. It is 
apparent that LR models do not adequately resolve the 
near-wall velocity gradient. However, different LR for- 
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mutations yield significantly different near-wall predic- 
tions as evidenced by the skin friction results. Of the four 
LR models examined, the SL technique yields the most 
accurate prediction for the backward-facing step flow in 
the 1:9 channel geometry. However, the standard HR k-c 
model is still superior due to the more accurate prediction 
of the skin friction behavior along the step side wall. 
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